Canonical Trajectories and Critical Coupling of the Bose-Hubbard Hamiltonian in a 

Harmonic Trap 
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Quantum Monte Carlo (QMC) simulations and the Local Density Approximation (LDA) are used 
to map the constant particle number (canonical) trajectories of the Bose Hubbard Hamiltonian 
confined in a harmonic trap onto the (/j,/U,t/U) phase diagram of the uniform system. Generically, 
these curves do not intercept the tips of the Mott insulator (MI) lobes of the uniform system. This 
observation necessitates a clarification of the appropriate comparison between critical couplings 
obtained in experiments on trapped systems with those obtained in QMC simulations. The density 
profiles and visibility are also obtained along these trajectories. Density profiles from QMC in the 
confined case are compared with LDA results. 
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The bosonic Hubbard model was first introduced^ in 
the context of disordered superconductors where the su- 
perfluidity of preformed Cooper pairs competes with 
Mott insulator and Bose glass phases. Considerable 
numerical work followed the original analytic treat- 
ment. When there is no disorder, Quantum Monte Carlo 
(QMC) studies^ obtained quantitative values for the 
critical coupling of the supcrfluid-Mott insulator (SF- 
MI) transition at commensurate filling in one dimen- 
sion, which were in good quantitative agreement with se- 
ries expansion^ and density matrix renormalization group 
calculations^. The critical point is now known in d = 2 
to a very high accuracy 7 -. 

Over the last decade, it became clear that trapped 
ultra-cold atoms provide an alternate, and more control- 
lable, experimental realization of the bosonic Hubbard 
model 8 -. Indeed, the possibility of a quantitative compar- 
ison of theoretical and experimental values for the critical 
point has been suggested. A recent experimental paper- 
has offered the first such benchmark in d = 2. 

However, a significant obstacle exists for such a direct 
comparison: The confining potential produces spatial in- 
homogeneities and a coexistence of SF and MI phases^. 
This naturally leads to the question as to what "criti- 
cal coupling" is being accessed in the experiments. Is it 
the coupling at which "Mott shoulders" begin to develop 
about a SF core? Or is it the coupling at which a Mott 
region pervades the entire central region of the trap? In 
this paper, we provide a detailed quantitative analysis 
of this issue. Specifically, using the Local Density Ap- 
proximation (LDA) and QMC simulations, we study, for 
fixed particle numbers, the evolution of the density pro- 
files of the trapped system as a function of the interaction 
strength and map those "canonical trajectories" onto the 
phase diagram of the uniform system. We also show data 
for the visibilit y 11 ' 12 . These measurements allow us to 
connect the critical points obtained in QMC with those 
that can be seen in experiment. 



The QMC results presented here were obtained using 
two different algorithms. In the firsts, the imaginary 
time /3 is discretized leading to a path integral for the par- 
tition function on a rigid space-imaginary time grid with 
local world line updates. In the secon d 16 ' 17 ' 18 , imaginary 
time is continuous and there are no Trotter errors associ- 
ated with discretization. Bosonic world-line updates can 
be non-local, and, as a consequence, the Green's function 
can be measured at all separations. The two algorithms 
give consistent results for all physical quantities calcu- 
lated such as the density profiles and superfluid density. 

The one dimensional bosonic Hubbard Hamiltonian is, 

H = -t (a\ai + i + o|+i°*) - 

i i 
i i 

Here i = 1, 2, ■ • • , L where L is the number of sites and 
Xi = a[i — L/2] is the coordinate of the ith site as mea- 
sured from the center of the system. We choose the lat- 
tice constant a = 1. The hopping parameter, t, sets 
the energy scale; in what follows we set t = 1, i.e., all 
energies are measured in units of t. m - 
number operator, and [a;,aj 
and destruction operators. Vt is the curvature of the 
trap, and the repulsive contact interaction is given by U . 
The chemical potential, fi, controls the average number 
of particles. 

The bosonic-Hubbard Hamiltonian can also be simu- 
lated in the canonical ensemble at fixed particle number 
Nf,. Indeed this is essential in order to make contact 
with experiments. In the homogeneous case, Vt = 0, 
the phase diagram is a function of the density, Nb/L d , 
and the interaction U jt where d is the dimensionality of 
the system. It was emphasized recently^ that a similar 
lattice size independent formulation can be made in the 
confined case by using a rescaled length £j = Xi/£ with 



a]a,i is the 



Sij are bosonic creation 
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£ = y/t/Vr- Then, density profiles and the resulting 
phase diagram depend on Nb and Vr only via the com- 
bination p = Nb/£ d , called the "characteristic density". 
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FIG. 1: (Color online) Density profiles i>s the rescaled posi- 
tion, in Id. The solid lines are obtained using QMC for 
the uniform system combined with the LDA to include the 
trap. The symbols are the results of QMC done directly on 
the confined system. The characteristic density p = 4.47 and 
U — 7.2,8.0,9.0. We also show, in the three panels, profiles 
for two different particle numbers, Nb = 50, 30 but with the 
same p = 4.47. 

One simple way to understand the role of the charac- 
teristic density and to infer the properties of the trapped 
system is the Local Density Approximation (LDA)^£ in 
which the density at a particular location Xi in the 
trapped system is assumed to be given by the density 
of a uniform system with chemical potential equal to the 
"local chemical potential" pi = p—Vrxf at that location. 



In other words, for a trapped Id system, 

p(xi) = {rii)v T = Pu((M> U), (2) 

where p° ld (p;U) = (rii)v T =o is the density for the Id 
bosonic Hubbard model in the homogeneous case. For 
a given desired Nb, the requisite chemical potential p in 
the presence of the trap, which is also the local chemical 
potential at the center of the trap, is determined by the 
condition, 



Nb = Jl(ni) VT =J2p°d(» - V T a$; U), 



(3) 



and is therefore implicitly a function of Nb, U and Vr- In 
principle, one can use p\ d {p] U) as determined by QMC 
in the uniform case, together with Eq. ([3]) to determine 
p. Within the LDA the density profile p(xi) is then com- 
pletely determined, and can be compared with results 
obtained directly from simulations with a trap potential 
to determine the accuracy of the LDA, as discussed be- 
low. Furthermore, Eq. ^ provides a useful guide to un- 
derstanding the trajectories in the (p, U) plane that are 
traversed in experimental investigations such as in Ref.i, 
since they are typically done at fixed Nb and varying t by 
varying the depth of the optical potential (which, how- 
ever, also changes the trap potential). 
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FIG. 2: (Color online) Canonical (constant particle number) 
flows in the (p/U,t/U) plane at fixed Vr = 0.008. Character- 
istic densities vary from p = 2.24 for Nb = 25 (lowest curve) 
to p = 10.73 for N b = 120 (highest curve). 

Better insight into the nature of such canonical ( con- 
stant Nb ) trajectories is obtained by approximating the 
sum in Eq. ([3]) as an integral. This should be a rea- 
sonable approximation when the local chemical potential 
changes slowly from site to site, i.e, in the same regime 
where LDA is expected to be valid. In the one dimen- 
sional case^ one has, 



N h = 2 



dxpl d (p-V T x 2 -U). 



(4) 



Changing the integration variable to p x = p — Vrx 2 , this 
equation can be rewritten as^, 



Nb\JV T = p ■ 



h{p;U). (5) 
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I\ is entirely determined from the solution of the homo- 
geneous problem. The chemical potential in the presence 
of the trap, p, is determined by inverting Eq. ([5|). Note 
the natural appearance of the characteristic density p on 
the left side of Eq. ([5]). Clearly, p, and hence the density 
profile expressed as a function of x/£, depend only on p 
and not on Nt, and Vt separately. Needless to say, I\ 
can also be computed directly by evaluating the sum in 
Eq. © using the simulation results for p® d . 

Thermodynamic stability implies that p\ d , and hence 
I\(p;U), are monotonically increasing functions^ of p. 
For p < Ii(ni(U);U) = p^(U), p, and hence pi, are less 
than p^(U), the chemical potential at which the first 
Mott lobe is reached from below. Therefore p(xi) < 1, 
and all sites are sampling the SF region in the phase 
diagram below the first Mott lobe (if U > U c ). 

The density profile is very different when U , Vt and Nf, 
are such that p is larger than p^(U). Then p > Pi(U) 
and therefore a flat Mott plateau with p(x{) = 1 appears 
in the central region of the system, extending over sites 
i for which pi > p^(U). For sites outside this plateau, 
p(xi) < 1 and the system is locally in the SF phase. 

If the trap potential is increased so as to squeeze the 
particles towards the center of the cell (or if JVf, is in- 
creased), p and p increase. For p > Ii(p^(U); U) = 
Pl(U), one has p > pf(U), the chemical potential at 
which the first Mott lobe is reached from above. In this 
case the central sites of the system are in the superfluid 
region above the Mott lobe, with p(xi) > 1, surrounded 
by MI shoulders where p(xi) = 1, in turn surrounded 
by SF regions as the edges of the system are reached 
(Fig. Ha)). 

In the regime p 1 (U) < p < p^(U), as is easily verified 
from Eq. ([5]), p is determined by the equation, 



p-p^(u) = [p- P r(u)] 2 /A. 



(6) 



Hence the two threshold values of p in the presence of the 
trap and the threshold chemical potentials for the Mott 
transition in the homogeneous case are related via, 



p+(U) - p^{U) = [p + (U) - j>r(U)¥/4 



(7) 



For larger values of p in large systems with a small 
Vt, one can access transitions involving the higher Mott 
lobes'. 

In Fig. [T] we compare the density profiles obtained from 
direct QMC simulations of the trapped system with those 
inferred from the LDA and QMC simulations of the uni- 
form system. The LDA generally provides an accurate 
description of the density profiles except at those loca- 
tions in the trap where a changeover from superfluid to 
Mott insulator region is occurring. This is clear in Fig. [1] 
where as one goes from SF to MI regions, the transition 
is much sharper for the LDA curves. This, of course, is 
a vestige of the true quantum phase transition present 
in the unconfined system on which the LDA method is 
based. Figure [1] shows profiles for two different pairs 
of (Nt,, Vt) which have the same characteristic density. 



They are seen to coincide almost perfectly, validating 
the use of and p to describe the physics in a scale- 
independent way. 

Figure [2] shows the canonical trajectories correspond- 
ing to p(p, U) for fixed p, obtained from Eq. ([3]), superim- 
posed on the phase diagram of the uniform system. Each 
trajectory is at constant Nf, and, therefore, constant p 
when Vt is fixed, and shows where the trapped system 
sits in the phase diagram of the uniform system when the 
LDA is used in combination with QMC. For example, for 
the confined system values Nb = 50, Vt — 0.008 and 
U = 9.0, p(p,U), lies well within the p = 1 Mott lobe 
and the system should be a Mott insulator according to 
this mapping. Figure [TJc) shows the true density profile 
obtained with QMC directly with a trap and we see that, 
indeed, the confined system is a MI, except for the edges 
which always have p[xi) < 1. On the other hand, staying 
on the same trajectory, Nb = 50, but with U = 7.2, the 
p(p, U) lies in the SF phase above the Mott lobe leading 
us to predict the central region of the trapped system to 
be SF with p(xi) > 1, as indeed confirmed by Fig. QJa). 
A second example of this evolution, for Nb = 110, which 
just clips the top of the p = 2 Mott lobe, is given in Fig.[3l 
with similar conclusions. Notice that as U increases, if a 
trajectory enters, say, the p = 2 Mott lobe, it will leave it 
eventually upon further increases in U. Such a trajectory 
will eventually enter the p — 1 Mott lobe which it can 
never leave. 

It is important to note that different trajectories in- 
tersect the Mott lobes at different (p/U,t/U) points and 
in general not at the tip. Thus, Fig. [2] emphasizes the 
central point of this paper, namely that both the parti- 
cle number and confining potential need to be considered 
together in determining the 'critical point' of the trapped 
boson Hubbard Hamiltonian. In particular, in order to 
access U c in an experiment, the characteristic density also 
has to be tuned to its appropriate critical value. In the 
case of a 1-d trapped system we are considering in this 
paper, p c ~ 2.7. 




FIG. 3: (Color online) Density profiles along the Nb = 110 
(p = 9.84) trajectory. This value just clips the tip of the 
uniform system p = 2 lobe, as seen in Fig. The dashed 
lines are to draw attention to the p = 1,2 values where the 
MI develops. 
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FIG. 4: (Color online) The visibility along the N b = 40 and 
Nb = 50 trajectories (corresponding to p = 3.58 and 4.47 with 
V T = 0.008). For N b = 50 the kink at U = 7.5 is associated 
with the presence of well-formed Mott shoulders. The second 
kink at U = 8.2 corresponds to the formation of a full Mott 
phase throughout the center of the trap (Vt — 0.008). 

Our understanding of the relation between the density 
profiles and the "flow diagram" of canonical trajectories 
is made complete by examining the visibility V, which is 
known to be a sensitive measure of the behavior of the 
density profiles 1 ^. For N b = 50 (p = 4.47), V has two 
kinks at U = 7.5 and U = 8.2 which indicate respectively 
the appearance of well-formed Mott shoulders surround- 
ing a SF interior and then the total disappearance of 
superfluidity at the trap center and the establishment of 
MI throughout (Fig. gj). It is seen from Fig. [2] that the 
second, larger, of these two values corresponds very well 
to the coupling where the 2V& = 50 trajectory enters the 



uniform system Mott lobe. 

In summary, in this paper we have shown that for fixed 
particle number, the "critical coupling" associated with 
destruction of superfluidity and onset of Mott behavior 
depends on the characteristic density p. In fact, this ob- 
servation is also implicit in the "state diagram" of [l(| 
in which the boundaries between phases at fixed Vt were 
shown to depend on JVj,. Using the local density approx- 
imation we explicitly constructed the trajectories in the 
(p/U, t/U) plane which correspond to constant p, and 
quantified their points of entry into the Mott lobe of the 
uniform system. This construction should allow experi- 
mentalists to predict where, on the phase diagram of the 
uniform system, their trapped system will be. The be- 
havior of the visibility confirmed that the uniform Mott 
lobe is entered when the center of the density profiles is 
in the Mott phase. 

We have focused here on d = 1. However, the basic 
qualitative point we wish to emphasize is valid in any 
dimension: a careful consideration of the confining po- 
tential in addition to the number of particles is essen- 
tial for a meaningful comparison of the critical couplings 
obtained in experiments with those of the homogeneous 
system. 
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